What is the effect of the teacher intervention on student grades among latinx students relative to their white peers?
Model levels:
LEVEL REMOVED: Data has been averaged over courses within teacher clusters (i.e., one row per student-by-teacher-by-cohort). Courses is considered a nuisance-level in this research context especially because course subject data is not available (e.g., science, math, English, etc.).
Imperfect hierarchy: The traditional Multi-Level Model (MLM) assumes that levels in the data are nested hierarchically. This dataset clearly violates this assumption given the students have multiple classes/teachers in a academic term. This data structure is described in quantitative methods literature as an imperfect hierarchy of the type multiple membership.
Two-way non-compliance: Although teachers were initially randomly assigned to the treatment conditions many teachers were unable to comply with their training assignment dates. Due to schedule conflicts or other unknown factors, teachers assigned to the treatment condition could not make the date of this training and then self-selected to the wait-list control condition. Teachers also self-selected into the treatment condition who were initially assigned to the control condition date. This situation is called two-way non-compliance in the causal inference literature. This issue occurs commonly in intervention studies, and although randomization is contaminated, it is still possible to recover a causal treatment estimate using a quasi-experimental method called instrumental variable analysis. In this case, the original assignment (assign) is the instrumental variable and the treatment condition received by a teacher is the binary variable treat.
We tried the IV modeling approach, but this method relies on a strong correlation between assignment and treatment. When the correlation between assignment and treatment is small, the instrument is said to be weak, which results in an invalid (biased) estimate of the causal effect. Therefore, other modeling approaches were considered.
CR/NC). Consequently a large proportion of the marks in this year were CR/NC which do not contribute to GPA and were coded as missing (NA).2019 COHORT DROPPED: This cohort was dropped based on the significant environmental events effecting grades in 2019. It was deemed that the grading policy changes resulted in grades with qualitatively distinct meaning relative to previous cohorts making grade comparisons across cohorts unfeasible.
mark1: Individual student-by-course grades (un-aggregated across teachers). Each student may have grades assigned by multiple teachers.
treat:
(0) Control = Teacher is in control group(1) Treatment = Teacher is in treatment groupwhite_latinx
(0) White = Student self-reported race is white(1) Latinx = Student self-reported race is latinxcohort: Year of intervention
STUDENT-LEVEL
school_type: Student’s school type is Junior High School (JH) or High School (HS)student_age: Students self-reported agesed: Student status is socio-economically disadvantaged (Not SED or SED)avg_absences: Average student absencesTEACHER-LEVEL
teach_white_latinx: Teacher is White, Latinx, or Othertitle: Teacher title is general classroom teacher (Teacher) or alternative teacher (Alt. Teacher)total_service_avg: Teacher’s total teaching years of serviceteacher_age: Teacher self-reported agePRIOR (REMOVED): Prior year average grades have been removed from analyses as a covariate. THis is due to challenges with access to this data as well as concerns with the prior year scores being contaminated by treatment.
## NULL
View all levels of course_subject
tabyl(droplevels(full_sample$course_subject))
library(DT)
datatable(full_sample, rownames = FALSE, filter="top",
options = list(pageLength = 15, scrollX=T) )
prop_treat)## full_sample$prop_treat n percent
## 0.0000000 4537 0.2214142794
## 0.1666667 42 0.0020496803
## 0.2000000 285 0.0139085452
## 0.2500000 1036 0.0505587819
## 0.2857143 28 0.0013664536
## 0.3333333 2232 0.1089258699
## 0.3750000 16 0.0007808306
## 0.4000000 640 0.0312332243
## 0.4285714 28 0.0013664536
## 0.5000000 4328 0.2112146796
## 0.5555556 9 0.0004392172
## 0.5714286 91 0.0044409741
## 0.6000000 725 0.0353813870
## 0.6250000 40 0.0019520765
## 0.6666667 1905 0.0929676443
## 0.7142857 168 0.0081987214
## 0.7500000 716 0.0349421697
## 0.8000000 425 0.0207408130
## 0.8333333 228 0.0111268362
## 0.8571429 77 0.0037577473
## 0.8750000 24 0.0011712459
## 1.0000000 2911 0.1420623688
tabyl(full_sample$course_type3)
## full_sample$course_type3 n percent
## non_acad 5212 0.2543556
## non_stem 7873 0.3842175
## stem 7406 0.3614270
Variable summary grouped by treatment status
full_sample %>%
dplyr::select(treat, cohort, mark1, course_type3,
teach_age, teach_white_latinx, teach_gender, school_type, number_students_16,
white_latinx, title, student_age, total_service_avg, sed, avg_absences) %>%
tbl_summary(by = treat,
statistic = list(all_continuous() ~ "{mean} ({sd})",
all_categorical() ~ "{n} ({p}%)")) %>%
add_n() %>% add_p() %>%
modify_header(label = "**Variable**") %>%
bold_labels()
| Variable | N | Control, N = 110061 | Treatment, N = 94851 | p-value2 |
|---|---|---|---|---|
| cohort | 20491 | <0.001 | ||
| 2017 | 7682 (70%) | 7709 (81%) | ||
| 2018 | 3324 (30%) | 1776 (19%) | ||
| mark1 | 20491 | 8.2 (3.4) | 8.5 (3.3) | <0.001 |
| course_type3 | 20491 | <0.001 | ||
| non_acad | 2167 (20%) | 3045 (32%) | ||
| non_stem | 3426 (31%) | 4447 (47%) | ||
| stem | 5413 (49%) | 1993 (21%) | ||
| teach_age | 20491 | 39 (11) | 43 (11) | <0.001 |
| teach_white_latinx | 20008 | <0.001 | ||
| White | 7579 (71%) | 7370 (79%) | ||
| Latinx | 1118 (10%) | 829 (8.9%) | ||
| Other | 2013 (19%) | 1099 (12%) | ||
| Unknown | 296 | 187 | ||
| teach_gender | 19993 | 0.066 | ||
| Male | 3568 (33%) | 3208 (35%) | ||
| Female | 7142 (67%) | 6075 (65%) | ||
| Unknown | 296 | 202 | ||
| school_type | 20491 | <0.001 | ||
| JH | 3881 (35%) | 3893 (41%) | ||
| HS | 7125 (65%) | 5592 (59%) | ||
| number_students_16 | 17039 | 26 (16) | 21 (11) | <0.001 |
| Unknown | 2354 | 1098 | ||
| white_latinx | 18709 | <0.001 | ||
| White | 4160 (41%) | 3971 (46%) | ||
| Latinx | 5920 (59%) | 4658 (54%) | ||
| Unknown | 926 | 856 | ||
| title | 20491 | <0.001 | ||
| Alt. Teacher | 517 (4.7%) | 206 (2.2%) | ||
| Teacher | 10489 (95%) | 9279 (98%) | ||
| student_age | 20489 | 13.92 (1.78) | 14.00 (1.82) | 0.008 |
| Unknown | 0 | 2 | ||
| total_service_avg | 20491 | 12 (9) | 16 (9) | <0.001 |
| sed | 20489 | <0.001 | ||
| Not SED | 5865 (53%) | 5508 (58%) | ||
| SED | 5141 (47%) | 3975 (42%) | ||
| Unknown | 0 | 2 | ||
| avg_absences | 20489 | 4.1 (4.5) | 4.0 (4.5) | 0.2 |
| Unknown | 0 | 2 | ||
|
1
Statistics presented: n (%); mean (SD)
2
Statistical tests performed: chi-square test of independence; Wilcoxon rank-sum test
|
||||
Multi-level modeling
library(psych)
library(rockchalk)
library(lme4)
library(varTestnlme)
library(lmerTest)
library(clubSandwich)
library(effects)
options(scipen = 999)
source(here("R-MLM-Functions.R"))
Teacher-level N = 65)treat: has variation at student & teacher level
white_latinx (level2): center within teacher?
If white_latinx is group-mean centered at the teacher-level could allow the between-level variable, within-level variable, or both to interact with treatment.
stem_prepared <- stem_sample %>%
mutate(white_latinx = as.numeric(white_latinx) - 1) %>%
mutate(cohort = as.numeric(cohort) - 1) %>%
mutate(school_type = as.numeric(school_type) - 1)
stem_centered <- gmc(stem_prepared,
c("white_latinx"),
by = c("teacher_id"),
FUN = mean,
suffix = c("_meanj", "_cwc"),
fulldataframe = TRUE) %>%
mutate(white_latinx_mj_cgm = white_latinx_meanj - mean(white_latinx_meanj, na.rm = TRUE)) %>%
mutate(cohort_cgm = cohort - mean(cohort, na.rm = TRUE)) %>%
mutate(school_type_cgm = school_type - mean(school_type, na.rm = TRUE))
Empty model with teacher-level & student-level random intercepts
mlm_level23 <- lmer(mark1 ~ 1 + (1 | teacher_id) + (1 | student_id), data = stem_sample, REML = FALSE)
mlm_level23
## Linear mixed model fit by maximum likelihood ['lmerModLmerTest']
## Formula: mark1 ~ 1 + (1 | teacher_id) + (1 | student_id)
## Data: stem_sample
## AIC BIC logLik deviance df.resid
## 37144.11 37171.75 -18568.06 37136.11 7402
## Random effects:
## Groups Name Std.Dev.
## student_id (Intercept) 2.313
## teacher_id (Intercept) 1.568
## Residual 2.050
## Number of obs: 7406, groups: student_id, 5330; teacher_id, 65
## Fixed Effects:
## (Intercept)
## 7.225
# ICC = (VAR_B1 + VAR_B2) / (VAR_B1 + VAR_B2 + VAR_W)
mlm <- as.data.frame(VarCorr(mlm_level23))
icc_students <- ((mlm[1,4] + mlm[2,4]) / (mlm[1,4] + mlm[2,4] + mlm[3,4]))
icc_teachers <- ((mlm[2,4]) / (mlm[1,4] + mlm[2,4] + mlm[3,4]))
# Proportion of variance captured by `student-level` random intercept:
round(icc_students*100,2)
## [1] 64.99
# Proportion of variance captured by `teacher-level` random intercept:
round(icc_teachers*100,2)
## [1] 20.46
white_latinx modelmlm_stem <- lmer(mark1 ~ 1 +
treat*white_latinx +
cohort_cgm +
school_type_cgm +
(1 | teacher_id) + (1 | student_id), data = stem_centered)
summary(mlm_stem)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mark1 ~ 1 + treat * white_latinx + cohort_cgm + school_type_cgm +
## (1 | teacher_id) + (1 | student_id)
## Data: stem_centered
##
## REML criterion at convergence: 33440.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.03092 -0.46349 0.06627 0.50531 2.77731
##
## Random effects:
## Groups Name Variance Std.Dev.
## student_id (Intercept) 4.942 2.223
## teacher_id (Intercept) 1.358 1.165
## Residual 4.059 2.015
## Number of obs: 6744, groups: student_id, 4874; teacher_id, 65
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 8.5516 0.1909 67.4460 44.796
## treatTreatment -0.4653 0.3503 65.7620 -1.328
## white_latinx -2.2317 0.0973 5752.1595 -22.937
## cohort_cgm -0.4566 0.3284 56.5169 -1.390
## school_type_cgm -1.2682 0.3382 56.1374 -3.749
## treatTreatment:white_latinx 0.5954 0.1635 4784.5016 3.641
## Pr(>|t|)
## (Intercept) < 0.0000000000000002 ***
## treatTreatment 0.188694
## white_latinx < 0.0000000000000002 ***
## cohort_cgm 0.169841
## school_type_cgm 0.000421 ***
## treatTreatment:white_latinx 0.000274 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtTrt wht_lt chrt_c schl__
## treatTrtmnt -0.526
## white_latnx -0.296 0.110
## cohort_cgm 0.015 -0.105 -0.010
## schl_typ_cg 0.020 -0.106 0.009 -0.203
## trtTrtmnt:_ 0.123 -0.279 -0.417 -0.011 -0.026
white_latinx within & between components interactmlm_stem <- lmer(mark1 ~ 1 +
treat*white_latinx_cwc +
treat*white_latinx_mj_cgm +
cohort_cgm +
school_type_cgm +
(1 | teacher_id) + (1 | student_id), data = stem_centered)
summary(mlm_stem)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mark1 ~ 1 + treat * white_latinx_cwc + treat * white_latinx_mj_cgm +
## cohort_cgm + school_type_cgm + (1 | teacher_id) + (1 | student_id)
## Data: stem_centered
##
## REML criterion at convergence: 33411.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.94604 -0.46685 0.06255 0.50382 2.76688
##
## Random effects:
## Groups Name Variance Std.Dev.
## student_id (Intercept) 4.8740 2.2077
## teacher_id (Intercept) 0.8274 0.9096
## Residual 4.1050 2.0261
## Number of obs: 6744, groups: student_id, 4874; teacher_id, 65
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 7.49198 0.14904 51.85344 50.267
## treatTreatment -0.08519 0.27628 49.82792 -0.308
## white_latinx_cwc -2.19194 0.09762 5717.32442 -22.454
## white_latinx_mj_cgm -5.50521 0.72585 54.22499 -7.585
## cohort_cgm -0.21044 0.26727 53.02994 -0.787
## school_type_cgm -1.25513 0.27743 52.79610 -4.524
## treatTreatment:white_latinx_cwc 0.62027 0.16494 4772.85664 3.761
## treatTreatment:white_latinx_mj_cgm 0.56014 1.17533 54.50462 0.477
## Pr(>|t|)
## (Intercept) < 0.0000000000000002 ***
## treatTreatment 0.759095
## white_latinx_cwc < 0.0000000000000002 ***
## white_latinx_mj_cgm 0.000000000452 ***
## cohort_cgm 0.434589
## school_type_cgm 0.000034739322 ***
## treatTreatment:white_latinx_cwc 0.000171 ***
## treatTreatment:white_latinx_mj_cgm 0.635569
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtTrt wht_l_ wht___ chrt_c schl__ trT:__
## treatTrtmnt -0.529
## wht_ltnx_cw -0.002 -0.003
## wht_ltnx_m_ -0.203 0.107 0.017
## cohort_cgm 0.041 -0.096 0.003 -0.124
## schl_typ_cg -0.001 -0.058 -0.002 0.138 -0.193
## trtTrtmn:__ 0.001 0.001 -0.416 -0.003 -0.001 -0.001
## trtTrtm:___ 0.116 -0.222 0.002 -0.615 -0.031 -0.223 0.002
Teacher-level N = 85)non_stem_prepared <- non_stem_sample %>%
mutate(white_latinx = as.numeric(white_latinx) - 1) %>%
mutate(cohort = as.numeric(cohort) - 1) %>%
mutate(school_type = as.numeric(school_type) - 1)
non_stem_centered <- gmc(non_stem_prepared,
c("white_latinx"),
by = c("teacher_id"),
FUN = mean,
suffix = c("_meanj", "_cwc"),
fulldataframe = TRUE) %>%
mutate(white_latinx_mj_cgm = white_latinx_meanj - mean(white_latinx_meanj, na.rm = TRUE)) %>%
mutate(cohort_cgm = cohort - mean(cohort, na.rm = TRUE)) %>%
mutate(school_type_cgm = school_type - mean(school_type, na.rm = TRUE))
Empty model with teacher-level & student-level random intercepts
mlm_level23 <- lmer(mark1 ~ 1 + (1 | teacher_id) + (1 | student_id), data = non_stem_sample, REML = FALSE)
mlm_level23
## Linear mixed model fit by maximum likelihood ['lmerModLmerTest']
## Formula: mark1 ~ 1 + (1 | teacher_id) + (1 | student_id)
## Data: non_stem_sample
## AIC BIC logLik deviance df.resid
## 38666.95 38694.83 -19329.47 38658.95 7869
## Random effects:
## Groups Name Std.Dev.
## student_id (Intercept) 2.416
## teacher_id (Intercept) 1.436
## Residual 1.772
## Number of obs: 7873, groups: student_id, 5775; teacher_id, 80
## Fixed Effects:
## (Intercept)
## 7.917
# ICC = (VAR_B1 + VAR_B2) / (VAR_B1 + VAR_B2 + VAR_W)
mlm <- as.data.frame(VarCorr(mlm_level23))
icc_students <- ((mlm[1,4] + mlm[2,4]) / (mlm[1,4] + mlm[2,4] + mlm[3,4]))
icc_teachers <- ((mlm[2,4]) / (mlm[1,4] + mlm[2,4] + mlm[3,4]))
# Proportion of variance captured by `student-level` random intercept:
round(icc_students*100,2)
## [1] 71.56
# Proportion of variance captured by `teacher-level` random intercept:
round(icc_teachers*100,2)
## [1] 18.67
white_latinx modelmlm_non_stem <- lmer(mark1 ~ 1 +
treat*white_latinx +
cohort_cgm +
school_type_cgm +
(1 | teacher_id) + (1 | student_id), data = non_stem_centered)
summary(mlm_non_stem)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mark1 ~ 1 + treat * white_latinx + cohort_cgm + school_type_cgm +
## (1 | teacher_id) + (1 | student_id)
## Data: non_stem_centered
##
## REML criterion at convergence: 34865.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8460 -0.3878 0.0840 0.4588 3.8396
##
## Random effects:
## Groups Name Variance Std.Dev.
## student_id (Intercept) 5.235 2.288
## teacher_id (Intercept) 1.622 1.273
## Residual 3.093 1.759
## Number of obs: 7188, groups: student_id, 5300; teacher_id, 80
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 9.17439 0.23919 94.31053 38.356
## treatTreatment -0.10874 0.31938 86.66579 -0.340
## white_latinx -2.01075 0.11646 7097.65153 -17.265
## cohort_cgm 0.01963 0.33671 73.58712 0.058
## school_type_cgm -0.97013 0.32118 75.24210 -3.021
## treatTreatment:white_latinx -0.07065 0.13851 6083.05070 -0.510
## Pr(>|t|)
## (Intercept) < 0.0000000000000002 ***
## treatTreatment 0.73432
## white_latinx < 0.0000000000000002 ***
## cohort_cgm 0.95368
## school_type_cgm 0.00345 **
## treatTreatment:white_latinx 0.61001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtTrt wht_lt chrt_c schl__
## treatTrtmnt -0.748
## white_latnx -0.312 0.195
## cohort_cgm -0.133 0.151 -0.004
## schl_typ_cg -0.168 0.174 0.004 -0.062
## trtTrtmnt:_ 0.223 -0.267 -0.704 -0.015 0.000
white_latinx (interaction for both within & between components)mlm_non_stem <- lmer(mark1 ~ 1 +
treat*white_latinx_cwc +
treat*white_latinx_mj_cgm +
cohort_cgm +
school_type_cgm +
(1 | teacher_id) + (1 | student_id), data = non_stem_centered)
summary(mlm_non_stem)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mark1 ~ 1 + treat * white_latinx_cwc + treat * white_latinx_mj_cgm +
## cohort_cgm + school_type_cgm + (1 | teacher_id) + (1 | student_id)
## Data: non_stem_centered
##
## REML criterion at convergence: 34853.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8338 -0.3886 0.0828 0.4595 3.8011
##
## Random effects:
## Groups Name Variance Std.Dev.
## student_id (Intercept) 5.221 2.285
## teacher_id (Intercept) 1.476 1.215
## Residual 3.101 1.761
## Number of obs: 7188, groups: student_id, 5300; teacher_id, 80
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 8.00568 0.22821 71.85462 35.081
## treatTreatment -0.13977 0.30274 70.89531 -0.462
## white_latinx_cwc -2.00775 0.11696 7044.93520 -17.166
## white_latinx_mj_cgm -2.24742 0.99173 75.95655 -2.266
## cohort_cgm 0.29689 0.33649 71.97193 0.882
## school_type_cgm -1.01251 0.30798 73.03861 -3.288
## treatTreatment:white_latinx_cwc -0.04858 0.13918 6013.27890 -0.349
## treatTreatment:white_latinx_mj_cgm -2.69726 1.40016 74.70703 -1.926
## Pr(>|t|)
## (Intercept) < 0.0000000000000002 ***
## treatTreatment 0.64573
## white_latinx_cwc < 0.0000000000000002 ***
## white_latinx_mj_cgm 0.02629 *
## cohort_cgm 0.38054
## school_type_cgm 0.00156 **
## treatTreatment:white_latinx_cwc 0.72709
## treatTreatment:white_latinx_mj_cgm 0.05786 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtTrt wht_l_ wht___ chrt_c schl__ trT:__
## treatTrtmnt -0.764
## wht_ltnx_cw -0.003 -0.002
## wht_ltnx_m_ -0.296 0.225 0.020
## cohort_cgm -0.153 0.162 -0.008 0.046
## schl_typ_cg -0.172 0.180 0.002 0.022 -0.071
## trtTrtmn:__ 0.002 0.002 -0.705 -0.006 0.007 -0.001
## trtTrtm:___ 0.234 -0.182 -0.002 -0.714 -0.232 0.014 0.000
Show NON-STEM & STEM model result side-by-side
library(sjPlot) # Latex MLM tables
tab_model(mlm_non_stem,mlm_stem,
dv.labels = c("NON-STEM","STEM"),
show.se = TRUE,
show.ci = FALSE
#p.val = "kr", # Kenward-Roger correction
#show.df = TRUE
)
| NON-STEM | STEM | |||||
|---|---|---|---|---|---|---|
| Predictors | Estimates | std. Error | p | Estimates | std. Error | p |
| (Intercept) | 8.01 | 0.23 | <0.001 | 7.49 | 0.15 | <0.001 |
| treat [Treatment] | -0.14 | 0.30 | 0.644 | -0.09 | 0.28 | 0.758 |
| white_latinx_cwc | -2.01 | 0.12 | <0.001 | -2.19 | 0.10 | <0.001 |
| white_latinx_mj_cgm | -2.25 | 0.99 | 0.023 | -5.51 | 0.73 | <0.001 |
| cohort_cgm | 0.30 | 0.34 | 0.378 | -0.21 | 0.27 | 0.431 |
| school_type_cgm | -1.01 | 0.31 | 0.001 | -1.26 | 0.28 | <0.001 |
|
treat [Treatment] * white_latinx_cwc |
-0.05 | 0.14 | 0.727 | 0.62 | 0.16 | <0.001 |
|
treat [Treatment] * white_latinx_mj_cgm |
-2.70 | 1.40 | 0.054 | 0.56 | 1.18 | 0.634 |
| Random Effects | ||||||
| σ2 | 3.10 | 4.11 | ||||
| τ00 | 5.22 student_id | 4.87 student_id | ||||
| 1.48 teacher_id | 0.83 teacher_id | |||||
| ICC | 0.68 | 0.58 | ||||
| N | 80 teacher_id | 65 teacher_id | ||||
| 5300 student_id | 4874 student_id | |||||
| Observations | 7188 | 6744 | ||||
| Marginal R2 / Conditional R2 | 0.144 / 0.729 | 0.202 / 0.666 | ||||